%% Carbon Taxes Around the World: Cooperation, Strategic Interactions, and Spillovers
% IMF Economic Review
% Alessandro Moro and Valerio Nispi Landi
% Replication files
% This file computes the initial steady state

function [ F ]  = find_steady(x, gamma,gammaz,eta,lam,bY,beta,betaz,GAM,GAMz,n,gY,gYz,FDIY,FDIYz,delta,iota,iY,iYz,...
                                 zeta,csi,mu,muz,nu,rkB,rkG,rkB_zu,rkB_zd,rkG_zu,rkG_zd,rkB_zz,chi,...
                                 rkG_zz,r,rz,price_clean,gdp_world,eWGtC,eGtC_share,Asym,gdp_sym,SYM,eme_factor_ab)
yH=x(1);
pH=x(2);
pB=x(3);
pBz=x(4);
zetae=x(5);
zetaez=x(6);
A=x(7);

pF=((1-(1-gamma)*pH^(1-eta))/gamma)^(1/(1-eta));
s=(gammaz*pH^(1-eta)+(1-gammaz)*pF^(1-eta))^(1/(1-eta));
gdp=pH*yH;
yFz=gdp/(lam*pF);
gdpz=pF*yFz/s;
b=bY*gdp;
dH=(1-beta/betaz)/GAM;
dF=s/GAMz*(betaz/beta-1);
v=-(b+dH+dF*(1-n)/n);
g=gY*yH;
gz=gYz*yFz;
k_zd=n*FDIY*gdp/(s*(1-n));
k_zu=FDIYz*gdpz*s*(1-n)/n;
i_zd=k_zd*(1-(1-delta)/iota);
i_zu=k_zu*(1-(1-delta)/iota);
i=iY*gdp-i_zu;
i_zz=iYz*gdpz-i_zd;
k=i/(1-(1-delta)/iota);
k_zz=i_zz/(1-(1-delta)/iota);
yB=zeta*(pB/pH)^(-csi)*yH;
conv_gdp=gdp_world/(n*gdp+(1-n)*pF*yFz);  
kappaA=price_clean*3.67/(conv_gdp*(n/zetae+(1-n)*s*eme_factor_ab/zetaez));
kappaAz=kappaA*eme_factor_ab;
tau=(mu^nu)*kappaA/zetae;                    
tauz=(muz^nu)*kappaAz/zetaez;                  
pB2=pB-tau*(1-mu)*zetae-kappaA/(1+nu)*mu^(1+nu);
pB2z=pBz-tauz*(1-muz)*zetaez-kappaAz/(1+nu)*muz^(1+nu);
pGY=pH*yH-pB*yB;
alpha1=k/(iota*(pB2/rkB*yB+pGY/rkG));
kB=alpha1*iota/rkB*pB2*yB;
kG=alpha1*iota/rkG*pGY;
alpha2=k_zu/(iota*(pB2/rkB_zu*yB+pGY/rkG_zu));
kB_zu=alpha2*iota/rkB_zu*pB2*yB;
kG_zu=alpha2*iota/rkG_zu*pGY;
hB=(yB/(A*(kB/iota)^alpha1*(kB_zu/iota)^alpha2))^(1/(1-alpha1-alpha2));
w=(1-alpha1-alpha2)*pB2*yB/hB;
hG=(1-alpha1-alpha2)*pGY/w;
h=hG+hB;
yG=A*(kG/iota)^alpha1*(kG_zu/iota)^alpha2*(hG)^(1-alpha1-alpha2);
pG=pGY/yG;
yBz=zeta*(s*pBz/pF)^(-csi)*yFz;
pGYz=gdpz-pBz*yBz;
alpha1z=k_zz/(iota*(pB2z/rkB_zz*yBz+pGYz/rkG_zz));
kB_zz=alpha1z*iota*pB2z*yBz/rkB_zz;
kG_zz=alpha1z*iota*pGYz/rkG_zz;
alpha2z=k_zd/(iota*(pB2z/rkB_zd*yBz+pGYz/rkG_zd));
kB_zd=alpha2z*iota*pB2z*yBz/rkB_zd;
kG_zd=alpha2z*iota*pGYz/rkG_zd;
c=-(i+pH*g+kappaA/(1+nu)*mu^(1+nu)*yB+b+(1-n)/n*s*i_zd-s*(1-n)/n*(rkB_zd*kB_zd    +rkG_zd*kG_zd)    /iota+(rkB_zu*kB_zu    +rkG_zu*kG_zu)    /iota+rz            *(dH+v)    /iota-r/iota    *(b+    dH    +v)-gdp);
cz=n/((1-n)*gammaz)*(pH/s)^(eta)*(yH-((1-gamma)*pH^(-eta)*(c+i+kappaA/(1+nu)*mu^(1+nu)*yB+i_zu)+g))-(i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz);
wz=(cz^(1/chi)*(1-alpha1z-alpha2z)*(pB2z*yBz+pGYz))^(chi/(1+chi));
hz=(wz/cz)^(1/chi);
hBz=(1-alpha1z-alpha2z)*pB2z*yBz/wz;
hGz=hz-hBz;
Az=yBz/((kB_zz/iota)^alpha1z*(kB_zd/iota)^alpha2z*(hBz)^(1-alpha1z-alpha2z));
yGz=Az*(kG_zz/iota)^alpha1z*(kG_zd/iota)^alpha2z*(hGz)^(1-alpha1z-alpha2z);
pGz=pGYz/yGz;
e=zetae*(1-mu)*yB;                                                                                                 
ez=zetaez*(1-muz)*yBz;   

if SYM==0

F=[h^(chi)*c-w
    (1-n)*((1-gammaz)*(pF/s)^(-eta)*(cz+i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz)+gz)+n*(gamma*(pF^(-eta))*(c+i+i_zu+kappaA/(1+nu)*mu^(1+nu)*yB))-(1-n)*yFz
    yGz-(1-zeta)*(s*pGz/pF)^(-csi)*yFz
    yG-(1-zeta)*(pG/pH)^(-csi)*yH
    n*e+(1-n)*ez-eWGtC
    n*e/eWGtC-eGtC_share
    n*gdp+(1-n)*gdpz-gdp_sym
    ];

elseif SYM==1
F=[h^(chi)*c-w
    (1-n)*((1-gammaz)*(pF/s)^(-eta)*(cz+i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz)+gz)+n*(gamma*(pF^(-eta))*(c+i+i_zu+kappaA/(1+nu)*mu^(1+nu)*yB))-(1-n)*yFz
    yGz-(1-zeta)*(s*pGz/pF)^(-csi)*yFz
    yG-(1-zeta)*(pG/pH)^(-csi)*yH
    n*e+(1-n)*ez-eWGtC
    n*e/eWGtC-eGtC_share
    A-Asym
    ];
end

end